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Abstract 

The two-dimensional t-J model in the ground state is investigated by the 
power Lanczos method. The pairing-pairing correlation function for d x 2_ y 2- 
wave symmetry is enhanced in the realistic parameter regime for high-T c 
superconductors. The charge susceptibility Xc shows divergent behavior as 
Xc oc near half-filling for the doping concentration 5, indicating that the 
value of the dynamical exponent z is four under the assumption of hyper- 
scaling. The peak height of the spin structure factor 5 max (Q) also behaves 
as S max (Q) oc near half-filling, which leads to the divergence of the an- 
tiferromagnetic correlation length £ m as £ m oc b^ 1 ! 2 . The boundary of phase 
separation is estimated on the basis of the Maxwell construction. Numerical 

results are compared with experimental features observed in high-T c cuprates. 
71.27.+a, 71.10.Fd, 71.30.+h, 75.40.Cx 
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I. INTRODUCTION 



Since the discovery of high-T c cuprate superconductors [Q], many microscopic models 
have been proposed in order to explain the pairing mechanism. The best way to under- 
stand the essential feature of the pairing mechanism would be to find the simplest and 
realistic model that describes the low-energy properties of the copper-oxide planes. The 
two-dimensional t-J model is one of the candidates for the effective model for the copper- 
oxide planes @||. As far as superconductivity is concerned, E. Dagotto and J. Riera have 
obtained indications of superconductivity by an unbiased diagonalization approach in the 
region of J/t ~ 3 near quarter- filling Q. More investigation is necessary in the realistic 
region of couplings and densities in order to confirm that the two-dimensional t-J model 
can really be an effective model for high-T c superconductors. One of the purposes of this 
paper is to investigate the relevancy of the two-dimensional t-J model as a low-temperature 
effective model for high-T c cuprates. 

Also, strong electron correlation in low-dimensional systems is one of the central issues 
in condensed matter physics. The Mott transition is one of the remarkable consequences 
of strong correlation. In general, there are two types of the Mott transitions for electron 
systems 0. One is characterized by the vanishment of the carrier density, the other is 
characterized by the divergence of the carrier mass. A typical example which shows the 
first type of transition is the free-fermion model on a lattice. In this case, the dynamical 
exponent z is two ||. A typical example which shows the second type of the Mott transition 
is the two-dimensional Hubbard model K has been shown numerically that the value 

of the dynamical exponent z is four in the case of the two-dimensional Hubbard model |||9| . 
A naive expectation is that the two-dimensional t-J model shows the same type of the Mott 
transition as that of the two-dimensional Hubbard model, because the t-J model can be 
derived as an effective model for the Hubbard model in the limit of U — > oo. However, in 
the large J/t regime, the t-J model shows different properties from those of the Hubbard 
model. For example, if J/t is so large that phase separation occurs, the transition to an 
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insulator is a first order transition. In the region of J/t slightly smaller than the phase- 
separation boundary, it is expected that electrons (or holes) form bound states due to the 
effective attractive forces that lead to phase separation. It is interesting to investigate the 
critical phenomena toward half-filling in this parameter regime. Hence, in this paper, we 
investigate the Mott transition of the two-dimensional t-J model in the ground state near 
half-filling. 

Several obstacles to numerical calculations have prevented us from getting the low-energy 
properties of the two-dimensional t-J model. For example, the system size achieved by the 
exact diagonalization is restricted to about 26 sites near half-filling. On the other hand, 
Quantum Monte Carlo algorithms have a serious sign problem. Recent progress of a Green's 
function Monte Carlo algorithm (power Lanczos method PI ) makes it possible for us to 
investigate the ground state properties of the t-J model in relatively large systems. In this 
paper, we use the usual Lanczos algorithm for clusters up to 20 sites and the power Lanczos 
method for larger clusters. 

In Sec. [Q], the t-J model is defined and the power Lanczos method is briefly reviewed. 
In Sec. [Til], we show the ground state energy as a function of filling and discuss phase 
separation on the basis of the Maxwell construction. In Sec. [TV], numerical results on 
the pairing-pairing correlation function are presented. In Sec. |V|, the Mott transition at 
J/t = 0.5 in the two-dimensional t-J model is discussed on the basis of the hyperscaling 



hypothesis. Section VI is devoted to summary. 



II. MODEL AND METHOD 

The t-J model is defined by the following Hamiltonian: 

7~CtJ = H t + Hj, 

'St 7, 



Ht = -t (Cqja + h.c.), 

<M><7 

<i,j> 



where cf CT denotes a creation operator of an electron at site i with spin a (a =T, !) with the 
constraint that no site is doubly occupied, which is defined as cf CT = (1 — ni^. a )cl a . The 
number operator n ia is defined as n ia = c\ a c i(T , using the standard electron creation operator 
cl a . The spin operator at site i is defined as Si = | Yl a p c L CT a/3Q/3, where cr^ is the vector 
of Pauli matrices. The summation (J2<i,j>) is taken over all nearest neighbor sites on a 
square lattice. 

We adopt the power Lanczos method proposed by Y.C. Chen and T.K. Lee ||10| . In the 
framework of this method, the expectation value of an operator O in the ground state of a 
Hamiltonian TC is evaluated by the following equation: 

(O) = lim(pLl|0|pLl)/ lim (pLl\pLl), (2.2) 

where |pLl) is the wavefunction defined as |pLl) = 7i p \Ll). The wavefunction | LI) is defined 
as |L1) = |trial) + Ci7Y|trial), where | trial) is a trial wavefunction and c% is a variational 
parameter, if we set c\ = 0, the power Lanczos method reduces to the power method. 



The variational wavefunction proposed by R. Valenti and C. Gros |IT] is employed as 
the trial wavefunction: 

10) = II l r i - rj-rPoP* |d-wave), (2.3) 

u 

where v is a variational parameter and represents the real-space coordinate at site i. The 
Gutzwiller and A^-particle projection operators are denoted by Pq and Pjv, respectively. 
The wavefunction |c?-wave) represents the BCS- wavefunction in which the order parameter 
has d x 2_ y 2-w&ve symmetry. It should be noted that the trial wavefunction is in the subspace 
that both the total spin S and the total momentum P are zero. Therefore the ground state 
properties reported in this paper are within this subspace. 

One of the reasons why we use the above wavefunction \<p) as a trial wavefunction is 
that this wavefunction gives the lowest energy among variational wavefunctions proposed so 
far, as far as we know, in the parameter regime we investigated. Another reason is that we 
can restrict the Hilbert space of simulation within the subspace of S = and P = 0. This 



makes convergence faster. The applicability of the power Lanczos method depends on the 
negative-sign ratio r which is defined by r = (p — n) / (p + n) , where p and n denote the 
number of positive and negative samples, respectively. If the ratio r is less than 0.1, it is 
difficult to obtain reliable results. As a result, the applicability of the power Lanczos method 
is restricted to the power p < p T , where p T denotes the power at which the ratio r is about 0.1. 
If the power required to reach convergence (p c ) is larger than p T , the power Lanczos method 
is not applicable. If we use a wavefunction which has small overlap with the ground-state 
wavefunction, it requires large power to reach the ground state. We compare the speed of 
convergence of the power Lanczos method and the simple power method using the Gutzwiller 
wavefunction and |0) as a trial wavefunction. As shown in Fig.[l], the power Lanczos method 
requires smaller power p than the simple power method, and the wavefunction |0) is superior 
to the Gutzwiller wavefunction as a trial wavefunction. In this figure, the error in energy 
due to finite power p is less than 0.2% for p > 5 by the power Lanczos method using 
although p c becomes larger than p T if we use the Gutzwiller wavefunction. 

We have checked convergent behavior in each simulation using |0) as a trial wavefunction. 
The p c becomes larger, if the system size becomes larger. For 50-site clusters, p c is about 
eight for the energy to converge. The p r becomes smaller, if J gets smaller. The most severe 
p x in our simulation is about eight near half-filling for J ~ 0.3. As an example, we show the 
convergence of energy in a 50-site cluster with 42 electrons at J = 0.3 in Fig.§. We measure 
physical quantities at p ~ 8, where we have checked in each simulation that the energy 
converges within a required accuracy. As a check of convergence of physical quantities, we 
show the pairing-pairing correlation function in a 20-site cluster with 18 electrons at p = 8 
in Fig.|| 

In the following sections, we show the numerical results in finite-size clusters up to 104 
sites. The boundary conditions are chosen for the momentum configuration to be closed- 
shell. We have typically run 1000-2000 Monte Carlo steps. Several hundred branches are 
produced at each Monte Carlo step in the evaluation of powers of 7i. 

The filling n is defined as n = N e /N B , where N e is the number of electrons and N s is 
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the number of sites. The doping concentration 5 is defined as 5 = 1 — n. The ground state 
energy per site at filling n is denoted by e(n). Hereafter we set t — 1 as the energy unit. 



III. PHASE SEPARATION 

Before investigating phase separation of the two-dimensional t-J model, we examine the 
finite-size effects on the ground state energy in the free-fermion model on a square lattice, in 
which we can calculate the exact ground state energy with any size of systems. We calculate 
the ground state energy per site of the free-fermion model in the same system sizes under the 
same boundary conditions as those used in the t-J model. We fit them as a function of filling 
by a polynomial up to third order. As shown in Fig.|], the fitting curve (dotted line) almost 
coincides with the curve in the thermodynamic limit (solid line). The finite-size effects on 
the ground state energy of the two-dimensional t-J model is probably not so different from 
those of the free-fermion model on a square lattice. Actually, as shown in Fig.|5], the data 
of the ground state energy per site of the two-dimensional t-J model in finite-size clusters 
are well fitted by a polynomial as a function of filling with small deviation from the fit, 
indicating that the finite-size effects on the ground state energy are small. 

Figure || shows the ground state energy per site as a function of filling in the two- 
dimensional t-J model at J = 0.5, 1.0, 1.5, 2.0, 2.5 and 3.0. In this figure, the tangent from 
the point at n — 1 to the fitting curve gives a lower energy than the fitting curve in the 
region of n c < n < 1 as represented by the solid line. Here n c is the electron density at 
the point of contact between the fitting curve and the tangent. Hence, we can identify the 
region of phase separation as n c < n < 1 on the basis of the Maxwell construction [12) . The 



energy and the chemical potential of the phase-separated state are given as the tangent and 
its slope, respectively. At J = 3.5, we find that the critical density n c is zero. On the other 
hand, as shown in Fig|| the chemical potential at J = 0.5 shows a monotonically increasing 
behavior as a function of filling at least in the region of n ~ 0.95, indicating that phase 
separation does not occur at J = 0.5. Hence, the phase-separation boundary is obtained as 
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in Fig.[12|. The critical J above which no stable homogeneous state exists is estimated as 
J C1 = 3.4 ± 0.1. The critical J below which no phase-separated state exists is estimated as 
J C2 = 0.75 ± 0.25. 

The critical value of J Cl is obtained more accurately in Ref.^ as 3.4367±0.0001 by solving 
the equation of motion of two electrons. The numerical result in this paper is consistent with 
it. The estimation of J C2 is consistent with that in Ref.^4|(J C2 ~ 0.75). In the intermediate 
region of n, the phase-separation boundary estimated in this paper is qualitatively similar 
to those in Ref.|] or Ref.|l5], but quantitatively lower than them (Fig.p^). 

In the following sections, we show the numerical results at J = 0.5. At this J, phase 
separation does not occur as discussed in this section. 



IV. PAIRING-PAIRING CORRELATION 

In this section, we show numerical results on the pairing-pairing correlation functions 
P±(r) defined as 

P ±(r) = W E( A ±(^o) t A ± (r + r)>. (4.1) 

iVs r-o 

Here, the singlet pairing operators A±(r) are defined as A(r) = c r j (c r+ £j + c r _f| ± c r+y i ± 
c r _yj), where + and — correspond to extended s-wave and <i x 2_ J/ 2-wave symmetry, respec- 
tively and the unit vectors in x- and y- directions are represented by x and y, respectively. 

In Fig. |6|(b), the pairing-pairing correlation function with d x 2_ y 2 -wave symmetry decays 
very little for n = 0.84, although the pairing-pairing correlation function with d x 2_ y 2-w&ve 
symmetry quickly decays for n = 0.20 as shown in Fig. |6](a). 

We define here the reduced pairing susceptibility as 

X± = E P±(r). (4.2) 

|r|>2 

Figures [7] (a) and (b) show the filling-dependence of x±/N s , where N s is defined as N s = 
J2\r\>2 1- If th e superconducting long-range order exists, the value of x±/N s remains finite in 
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the thermodynamic limit. In Figs.|7| (a) and (b), the pairing-pairing correlation for d x 2_ y 2- 
wave symmetry is enhanced in the region of 0.6 ~ n ~ 1, and that for extended s-wave 
symmetry is a little enhanced in the low-density regime. 

The numerical results showing that the d^.^-wave component of the pairing-pairing 
correlation is dominant near half-filling are consistent with experimental indications, for 
example, the measurements of the phase coherence in bimetallic YBCO-Pb dc SQUIDs |L6 



V. MOTT TRANSITION 

In Fig. §, the filling- dependence for the chemical potential at J = 0.5 is shown. The 
data of the chemical potential in finite-size clusters are calculated as follows: 

_ e(wi) -e(w 2 ) ( . 

rii - n 2 

where n is taken as n = (rii + n 2 )/2. Here, n\ and n 2 are taken to be adjacent closed-shell 
filling with boundary conditions fixed. Here, boundary conditions at half-filling are regarded 
as those under which the momentum configurations are closed-shell in the free-fermion model 
on a square lattice. In the thermodynamic limit, this definition of the chemical potential 
reduces to the normal one: fi(n) = de(n)/dn. 

We fit the data near half-filling in Fig||(c) as fi — fi c oc 5 a and estimate /i c = 1.31 ± 0.03 
and a = 1.78 ± 0.29, which is close to a = 2 reported in the case of the two-dimensional 
Hubbard model at U = 40. This suggests that the charge susceptibility Xc defined by 
Xc = dn/dfi diverges as Xc oc S^ 1 toward half-filling. Actually the chemical potential is 
fitted well by the following form: 

|/i — fx c \ oc 5 2 , (5.2) 

as denoted by the dashed line in Fig.|8](a) and (b). In order to check this divergent behavior 
of the charge susceptibility, we also investigate the doping dependence of chemical potential 
for J = 0.3 and J = 0.4. Figure |5] shows the same plot as in Fig.|S] for J = 0.3 and J = 0.4. 
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From the fit of the numerical data in Fig||(c), we estimate /x c = 1.88 ±0.03, a = 1.83 ±0.17 
for J = 0.3 and /i c = 1.57 ± 0.02, a = 1.98 ± 0.07 for J = 0.4. The numerical results for 
J = 0.3 and J = 0.4 also suggest that the charge susceptibility Xc diverges as Xc oc <5 _1 . 
The divergent behavior of the charge susceptibility is consistent with recent photo-emission 
measurements |I7| . 



If the hyperscaling relations are satisfied, the charge susceptibility Xc near the transition 
point to an insulator is written as Xc oc where z is the dynamical exponent and d 

is the spatial dimensionality Q. Therefore the numerical results suggest that the value of 
the dynamical exponent z is four. 

We define the spin structure factor S(k) as 

S(k) = W( S o'S r )e ikr . (5.3) 

6 r 

Figure |10|(a) shows the peak height of the spin structure factor 5' max (Q) as a function of 
filling. The data near half-filling can be fitted well by the following form: 

s m UQ) cx r 1 , (5.4) 

as denoted by the dashed line(Fig.|TD|(a), (b) and (c)). We fit the data near half- filling as 
5'max(<5) _1 and estimate (3 = 1.02 ± 0.02. This suggests that the antiferromagnetic 

correlation length £ m diverges toward half-filling as 

U oc r 1 / 2 , (5.5) 

under the assumption that the spin-spin correlation behaves as (So • S r ) oc e l( ^ r ■ e~ r ^ m 0. 
This behavior of the correlation length has been reported on the two-dimensional Hubbard 
model at U — 4 |7[] and is consistent with the observation by neutron scattering experiments 

Under the assumption of the existence of the single characteristic length scale £ that is 
related to critical phenomena, the hyperscaling theory has predicted that the length scale 
£ diverges as £ oc S^ 1 ^' toward the critical point, where d is the spatial dimensionality |§. 
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The numerical results shown above support the scaling hypothesis and suggest that the Mott 
transition in the two-dimensional t-J model at J = 0.5 is characterized by the dynamical 
exponent 2 = 4, which is the same as in the case of the two-dimensional Hubbard model at 

VI. SUMMARY 

Numerical results presented in this paper are consistent with the following experimental 
features found in the high-T c oxides: (i) the (i x 2_ J/ 2-wave symmetry of the superconducting 
order parameter in the region of moderate doping, which is suggested by the measurements 
of the phase coherence in bimetallic YBCO-Pb dc SQUIDs |Tj| (Sec. |TV| ), (ii) the doping 
dependence of the antiferromagnetic correlation length near half-filling (£ m oc S^ 1 ^ 2 ) ob- 
served in neutron scattering experiments [|18| (Sec. |V|), (iii) the large Fermi surface behavior 
in the region of moderate doping ]T{5| (Fig.|TT|) and (iv) divergent behavior of the charge 
susceptibility suggested by photo-emission experiments |L7| (Sec. |V|). 

In summary, numerical results on the two-dimensional t-J model have been reported. 
The boundary of phase separation is estimated on the basis of the Maxwell construction 
(Fig.|l!g). The pairing-pairing correlation for dj.2_ y 2-wave symmetry is enhanced in the region 
of 0.6 ~ n ~ 1 at J = 0.5 (Fig.[7|). The charge susceptibility Xc shows divergent behavior 
as Xc oc S^ 1 toward half-filling, indicating that the value of the dynamical exponent z is 
four (Figs.|| and |9|). The peak height of the spin structure factor S max (Q) diverges toward 
half-filling as S max ,(Q) oc <5 _1 (FigjT0|). This leads to the divergence of the antiferromagnetic 
correlation length as £ m oc b~ x l 2 . 
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FIGURES 

FIG. 1. Energy per site as a function of power p for the 20-site system with 18 electrons at 
J = 0.5. The percentage of the error from the ground state energy corresponds to the vertical scale 
on the right. Open circles and open diamonds denote the data obtained by the power method and 
the power Lanczos method, respectively, using the Gutzwiller wavefunction as a trial wavefunction. 
Solid symbols denote the data by using \(j>) as a trial wavefunction. 

FIG. 2. The same plot as in Fig.[l] but for the 50-site system with 42 electrons at J = 0.3. Inset 
shows the negative-sign ratio r defined in the text. The dashed line represents r = 0.1. 

FIG. 3. Pairing-pairing correlation function for d^^-wave symmetry in a 20-site cluster with 
18 electrons at J = 0.5. Crosses denote the data by \4>) without the power method. Solid diamonds 
denote the data by the power Lanczos method at p = 8, using \cj>) as a trial wavefunction. Open 
circles denote the exact data obtained by the exact diagonalization. 

FIG. 4. Ground state energy per site as a function of filling in the two-dimensional free-fermion 
model on a square lattice. Dotted line denotes a polynomial fit. Solid line denotes the ground state 
energy per site in the thermodynamic limit. 

FIG. 5. Ground state energy per site as a function of filling in the two-dimensional t-J model 
at J = 0.5,1.0,1.5,2.0,2.5 and 3.0, starting from above. Dotted line denotes the same fit as in 
FigJ^, using data points from n = to n ~ 0.7. Solid line denotes the expected ground state energy 
per site of the phase-separated state which is determined on the basis of the Maxwell construction. 
Dashed line at J = 0.5 is obtained by integrating the fit (dashed line) in Fig.||. 

FIG. 6. Pairing-pairing correlation function as a function of distance for (a) n = 0.20 and (b) 
n = 0.84 at J = 0.5 in 50-site clusters. Crosses and solid diamonds denote the pairing-pairing 
correlation functions for extended s-wave symmetry and d x 2_ y 2-w&ve symmetry, respectively. 
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FIG. 7. Reduced pairing susceptibility per site (x±/N s ) as a function of filling at J = 0.5 
for (a) extended s-wave symmetry and (b) d x 2_ y 2-wave symmetry, where x± = J2\r\>2 P±( r ) an d 
^s = E|r|> 2 l- 

FIG. 8. Chemical potential as a function of filling at J = 0.5, (a) linear plot, (b) fi vs S 2 plot 
and (c) log-log plot. Dotted line in (a) and (b) is obtained by differentiating the fit at J = 0.5 
(dotted line) in Fig.||. Dashed line in (a) and (b) denotes a fit as |/i — fi c \ oc 5 2 . Dashed and solid 
lines in (c) correspond to the cases of the dynamical exponent z = 2 and z = 4, respectively. 

FIG. 9. Chemical potential as a function of filling at J = 0.3, 0.4 and 0.5 starting from above, 
(a) linear plot, (b) /J, vs 5 2 plot. (c)Log-log plot for J = 0.3 (open symbols) and J = 0.4 (solid 
symbols). Dashed line in (a) and (b) denotes a fit as \{i — fi c \ oc 5 2 . The point denoted by a cross 
is obtained by a Green's function Monte Carlo in a lOxlO-site cluster with two holes taken from 
Ref.|2^. Dashed and solid lines in (c) correspond to the cases of the dynamical exponent z = 2 and 
z = 4, respectively. 

FIG. 10. Peak height of the spin structure factor S max (Q) as a function of filling at J = 0.5, 
(a) linear plot, (b) 5 max (Q) _1 vs 5 plot and (c) log-log plot. Dashed line in (a), (b) and (c) denotes 
a fit in the form of S max (Q) oc 

FIG. 11. Momentum distribution function for n = 0.84 at J = 0.5 in a 50-site cluster. In the 
inset, solid and dashed lines denote the Fermi Surface and the Brillouin zone boundary, respectively. 
The center is the T-point. 

FIG. 12. Schematic phase diagram of the two-dimensional t-J model in the ground state. Solid 



diamonds denote the phase-separation boundary obtained in Sec. III. In the higher-density region 
than the open diamond, the pairing-paring correlation for ci a .2_ 2/ 2-wave symmetry is enhanced. It 
should be noted that a ferromagnetic phase may exist in the small J region (J ~ 0.1) suggested in 
Refs.|l5| and 21. This ferromagnetic phase is beyond the scope of this paper. 
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